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We investigate the phenomenon of extinction of a long-lived self-regulating stochastic population, 
caused by intrinsic (demographic) noise. Extinction typically occurs via one of two scenarios de- 
pending on whether the absorbing state n = is a repelling (scenario A) or attracting (scenario B) 
point of the deterministic rate equation. In scenario A the metastable stochastic population resides 
in the vicinity of an attracting fixed point next to the repelling point n = 0. In scenario B there is an 
intermediate repelling point n = Ui between the attracting point n — and another attracting point 
n = U2 in the vicinity of which the metastable population resides. The crux of the theory is a dissi- 
pative variant of WKB (Wentzel-Kramers-Brillouin) approximation which assumes that the typical 
population size in the metastable state is large. Starting from the master equation, we calculate 
the quasi-stationary probability distribution of the population sizes and the (exponentially long) 
mean time to extinction for each of the two scenarios. When necessary, the WKB approximation 
is complemented (i) by a recursive solution of the quasi-stationary master equation at small n and 
(ii) by the van Kampen system-size expansion, valid near the fixed points of the deterministic rate 
equation. The theory yields both entropic barriers to extinction and pre-exponential factors, and 
holds for a general set of multi-step processes when detailed balance is broken. The results simplify 
considerably for single-step processes and near the characteristic bifurcations of scenarios A and B. 

PACS numbers: 02.50.Ga, 87.23. Cc 



I. INTRODUCTION 

Extinction of an isolated stochastic population after 
maintaining a long-lived state is a dramatic phenomenon. 
It occurs, even in the absence of environmental varia- 
tions, because of an unusual chain of random events when 
population losses dominate over gains. Population ex- 
tinction risk is a key negative factor in viability of small 
populations [l|, 0] , whereas extinction of a disease follow- 
ing an epidemic outburst is of course favorable. The 
possibility and consequences of extinction of biologically 
important components, regulated by chemical reactions 
in living cells, have also attracted interest As stochas- 
tic population dynamics are usually far from equilibrium, 
and no general methods of evaluating large fluctuations 
are available, they are of much interest to physics HQ. 

This work deals with an isolated single-species popu- 
lation undergoing a set of gain-loss processes. We will 
assume that the population is well mixed, so that spa- 
tial degrees of freedom are irrelevant. At the level of the 
deterministic rate equation (henceforth rate equation), 
which describes the time history of the mean population 
size n(t) and ignores fluctuations, h(t) flows to an attract- 
ing fixed point, where the gain and loss processes balance 
each other. The actual stochastic population, however, 
behaves differently and ultimately becomes extinct. This 
is because, in the absence of influx of new individuals, the 
empty state n = is absorbing: the probability of exiting 
from it is zero [7[. 

Although extinction (and fluctuations in general) are 
beyond its scope, the rate equation is a convenient start- 
ing point of our analysis. For an isolated single-species 
population the rate equation can be written as 

dn 

— =n#(n), (1) 
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FIG. 1: (color online). Typical extinction scenarios are de- 
termined by the character of the fixed point n = of the rate 
equation |T]). (a) Scenario A: the fixed point n — is repelling. 
In the stochastic system extinction occurs via a large fluctua- 
tion which brings the metastable population from a vicinity of 
the attracting fixed point n = m of the rate equation directly 
to the absorbing state n = 0. (b) Scenario B: the fixed point 
n — is attracting. In the stochastic system extinction occurs 
via a large fluctuation which brings the metastable popula- 
tion from a vicinity of the next attracting fixed point n = 712 
of the rate equation to a vicinity of the repelling fixed point 
n = in. From there the population flows "downhill" to the 
absorbing state n = almost deterministically. 

where is a smooth function determined by the spe- 
cific gain-loss processes, see below. For generic gain-loss 
processes $'(0) 7^ 0. For $'(0) > the fixed point fi = 
is repelling, whereas for $'(0) < it is attracting. In the 
former case, the next fixed point h = n\ > of Eq. (p} is 
attracting, see Fig. [Ij,. According to the rate equation, 
the mean population size in this case flows to n = n\ and 
stays there forever. When varying the rate constants of 
the gain-loss processes, the attracting fixed point n = n\ 
emerges via a transcritical bifurcation. 

Now let n = be an attracting fixed point of the 
rate equation ((T|). To have a long-lived population of 
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a nonzero size, at least two more fixed points of the 
rate equation (JTJ) must be present: a repelling point 
n = tlx > and an attracting point n — n 2 > n±, see 
Fig. [TJd. When starting from any n[t = 0) > m, the 
mean population size flows to n — n 2 and, according to 
the rate equation, stays there forever. The characteristic 
bifurcation in this case is saddle-node. 

As we will see shortly, these two cases give rise to two 
different extinction scenarios of stochastic populations. 
To account for the intrinsic noise, we employ the master 
equation 

= [W r (n - r)P n _ r (t) - W r (n)P n (t)] (2) 

r 

which describes the evolution of the probability P n (t) to 
have n individuals at time t. Here W r (k) > is the 
transition rate between the states with k and k + r in- 
dividuals, whereas r = ±f,±2,..., and all terms that 
include Pk with k < are assumed to be zero. For Po(t) 
the master equation is 



' r<0 



(3) 



For n = to be an absorbing state, the process rates must 
obey, for any r = ±1, ±2, . . . , the condition Wi-(0) = 0. 

We will be interested in the important regime of pa- 
rameters for which the mean population size in the 
metastable state, as predicted by Eq. ((T|), is large com- 
pared to one. Here, prior to extinction, a long-lived prob- 
ability distribution function (PDF) of the population sets 
in, on a relaxation time scale t r , around the correspond- 
ing attracting fixed point of the rate equation. This long- 
lived PDF, however, is metastable: it slowly decays in 
time. Simultaneously, the probability to find the popu- 
lation extinct slowly grows in time, see e.g. Refs. [8|, [9(: 



Pn>o{t > *r) ^ 7T„e * /t , P (t > t r ) ~ 1 - e 



t/7 



(4) 



The shape function 7r„ (n = 1,2,...) of the metastable 
PDF is called the quasi-stationary distribution (QSD). 
For metastable populations a very strong inequality, 
t ^> t r holds, and the decay time r is equal to the 
mean time to extinction (MTE): the mean time it takes 
the stochastic process to reach the absorbing state at 
n = 0. The main objectives of this work is to accurately, 
and analytically, calculate the QSD ir n and the MTE r 
of a population which experiences quite a general set of 
stochastic gain-loss processes. The crux of the method 
is a dissipative WKB approximation fl(il-[l2l]. where one 
assumes n ^> 1, treats n as a continuous variable and 
searches for 7r„ as 



-NS(n)-S 1 (n)-(l/N)S 2 (n)- 



(5) 



Here N 1 is a large parameter which scales as the 
mean population size in the metastable state. S(n) is 
called the action, whereas a(n) = e~ Sl (") is called the 



FIG. 2: (color online). Example of scenario A of popula- 
tion extinction driven by intrinsic noise. Shown are the zero- 
energy trajectories of the WKB Hamiltonian H{n,p) for the 

reactions A A 2 A, A A- and 2 A A 16|. The trajectories 
denoted by the thicker line determine the WKB solution for 
7r n , obtained in Ref. [l6|. The activation trajectory connects 
the metastable point (ni , 0) and the fluctuational extinction 
point (0,pf), where pf = — lni?o, and Ro = A/^i. The ef- 
fective entropy barrier to extinction is equal to NAS, where 
N — X/a and AS is the area of the shaded region, given by 
Eq. (7BJ|. 



amplitude. The WKB approximation breaks down at 
n = 0(1). Here a different approximation must be used, 
as explained below. 

Here is an overview of the two extinction scenarios as 
described by the WKB approximation. First, let n — 
be a repelling fixed point of the rate equation, see Fig.QJi. 
In a stochastic description extinction occurs via a large 
fluctuation which, acting against an effective entropy bar- 
rier, brings the population from a vicinity of n = n\ 
directly to the absorbing state n = 0. In the WKB lan- 
guage this transition is possible because of the presence of 
the fluctuational momentum p = dS/dn, see Fig. [21 The 
attracting and repelling fixed points of the rate equa- 
tion n = n\ and n = 0, respectively, become hyperbolic 
fixed points of an extended phase plane (n,p). Impor- 
tantly, an additional hyperbolic fixed point (0,pf) - the 
fluctuational extinction point - appears here, with a zero 
coordinate, n = 0, but a nonzero momentum pf 
IToT ]. The most probable path to extinction is the het- 
eroclinic trajectory, directly connecting the "metastable 
point" , that is the hyperbolic point (ni, 0), and the "fluc- 
tuational extinction point": the hyperbolic point (0,pf). 
(Such escape trajectories - heteroclinic trajectories with 
a non-zero momentum - are often called "activation tra- 
jectories", see e.g. [H].) This is what we call extinction 
scenario A. 

Now let n = be an attracting fixed point of the rate 
equation (£Q), so that the metastable population resides in 
the vicinity of n = ri2, see Fig. [TJd. In the stochastic de- 
scription, extinction occurs via a large fluctuation which 
brings the population from a vicinity of n — n% to a vicin- 
ity of the repelling fixed point n± . From there the system 
flows into the absorbing state n = "downhill", that is 
almost deterministically. In the framework of WKB the- 
ory the transition from n 2 to n\ occurs in the extended 
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FIG. 3: (color online). Example of scenario B of popula- 
tion extinction driven by intrinsic noise. Shown are the zero- 
energy trajectories (|106[) on the phase plane (n, p) of the WKB 

Hamiltonian (|105[) for the reactions A A 0, 2 A A 3 A and 
3 A A 2 A. The trajectories denoted by the thicker lines de- 
termine the WKB solution for the QSD. The most probable 
path to extinction first goes along the non-zero-momentum 
heteroclinic trajectory (the activation trajectory) which con- 
nects the points (712, 0) and (ni, 0). Then the population flows 
almost deterministically to the extinction point (0, 0) along a 
zero-momentum segment (the relaxation trajectory). The ef- 
fective entropy barrier to extinction is equal to NAS, where 
N = 3A/(2a), and AS is the area of the shaded region, given 
by Eq. (fTlO)) . 



phase plane (n, p) where all three fixed points are hyper- 
bolic, see Fig. [3] Here the optimal path to extinction 
is composed of two segments: the non-zero-momentum 
heteroclinic trajectory connecting the hyperbolic fixed 
points (^2,0) and (nj.,0) (the activation trajectory), and 
the zero-momentum segment going from n — ri\ to n = 
(the relaxation trajectory). This is what we call extinc- 
tion scenario B. 

The mean time to extinction (MTE) r and/or the 
QSD ir n of metastable single-species stochastic popula- 
tions were calculated previously in particular examples 
in different contexts of physics, chemistry, population bi- 
ology, epidemiology, cell biology, etc. Among them there 
is a large body of work which approximated the master 
equation by an effective Fokker-Planck equation, derived 
via the van Kampen system size expansion or related 
recipes. Once the Fokker-Planck equation is obtained, 
the MTE and QSD can be calculated by standard meth- 
ods Unfortunately, this approximation is in general 
uncontrolled. It fails in its description of the tails of the 
QSD, and gives exponentially large errors in the MTE, 
as shown in Refs. (ig-[l9j. 

With a few exceptions, accurate analytic results for the 
MTE and QSD are only available for single-step gain- 
loss processes: r = ±1 in Eq. ©. In this case the MTE 
can be determined exactly by employing the backward 
master equation [5|, |6|. This yields a cumbersome an- 
alytic expression for the MTE which, for a large pop- 
ulation size in the metastable state, can be simplified 
via a saddle-point approximation. Such a procedure was 
implemented in Ref. [l8| . In its turn, the QSD tt„ of 
single-step processes can be calculated from a recursive 



relation obtained when substituting Eq. Q in the master 
equation. Several model examples of multi-step processes 
were considered in Refs. d, [T(| El) H3|, all of them be- 
longing to extinction scenario A. We will generalize the 
previous results substantially and determine the MTE 
and QSD for quite a general set of gain-loss processes 
pertaining to extinction scenario A. We will also deter- 
mine the MTE and QSD for extinction scenario B. 

Our WKB theory starts with applying the ansatz ([5]) 
to an eigenvalue problem for the QSD 7r„ which is nothing 
but the first excited eigenvector of the master equation. 
In the leading WKB order one arrives at the problem 
of finding zero-energy trajectories of an effective classical 
Hamiltonian [121 ] . There are two different types of zero- 
energy phase trajectories (in addition to the extinction 
line q = 0): the activation and relaxation trajectories, 
which correspond to the fast and slow WKB modes, re- 
spectively [21] . To obtain the pre-exponents, one needs 
to consider the sub-leading WKB order. The WKB cal- 
culations are simpler for scenario A, as the relaxation 
trajectory does not play any role here. In scenario B 
both the activation, and the relaxation trajectories are 
important. In both scenarios the WKB approximation 
breaks down at n = 0(1). Here we find the QSD, up to 
a normalization constant, from a recursive relation, ob- 
tained by linearizing the process rates with respect to n 
at sufficiently small n. In scenario A it suffices to match 
the recursive solution with the fast-mode solution in their 
joint region of validity, in much the same way as it was 
done by Kessler and Shnerb [l6[ in a particular exam- 
ple of three stochastic reactions. In Scenario B the slow 
mode dominates the WKB-solution at n < n\. It di- 
verges, however, at n = n\. To obtain a regular solution 
there, one needs to go beyond the WKB approximation 
and account, in a close vicinity of n — n\, for strong cou- 
pling between the fast- and slow-mode solutions. This 
can be done via the van-Kampen system size expansion 
of the master equation which does hold in the vicinity 
of n = m. This procedure was first implemented by 
Meerson and Sasorov [2l[ , in a model problem of noise- 
driven population explosion. Then it has been employed 
by Escudero and Kamenev [22| in the context of a WKB 
theory of stochastic population switches between two dif- 
ferent metastable states. The theory of Escudero and 
Kamenev [22j was formulated for quite a general set of 
gain-loss processes. In this paper we will adopt their gen- 
eral approach, and some of their notation, in the problem 
of population extinction. 

Here is a plan of the remainder of the paper. Section HT1 
starts with a formulation of the eigenvalue problem for 
the QSD. Then we expose the WKB approximation and 
the fast- and slow-mode WKB solutions. The derivation 
here is quite general and holds for extinction scenarios 
A and B. Section IIIII presents a derivation of recursive 
solution of the quasi-stationary master equation for suf- 
ficiently small n. This derivation, which also holds for 
extinction scenarios A and B, is specific to population ex- 
tinction. Except for simple particular cases (see e.g. Ref. 
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[IB]), it has not been attempted before. In Section HVl we 
match the fast-mode WKB solution with the recursive 
small-n solution and obtain general expressions for the 
QSD and MTE in scenario A. In the same Section we ob- 
tain the QSD and MTE for single-step processes and near 
the transcritical bifurcation, characteristic of scenario A. 
Then we illustrate our theory on several particular exam- 
ples, some of which investigated previously. In Section [V] 
we determine the QSD and MTE for scenario B. Then we 
again apply our results to single-step processes and near 
the saddle- node bifurcation, characteristic of scenario B. 
Furthermore, we consider a particular example of three 
stochastic reactions and compare our theoretical predic- 
tions with a numerical solution of the master equation. 
A summary of our results is presented in Section IVI1 



and is crucial both for the WKB-approximation that we 
present in this Section, and for the recursive solution of 
Eq. ([7]) that we will be dealing with later. As q = is 
the absorbing state, w r (0) = u r (0) — 0. 

For n>lwe can employ the WKB ansatz ([5]): 



n(q) = ir Nq 



Ae 



-NS(q)-S 1 (q) 



(10) 



where S(q) and S±(q) are assumed to be 0(1), and a 
constant prefactor A is introduced for convenience, see 
below. Now we assume that |r mQ2; | -C N, Taylor-expand 
the functions of q — r/N in Eq. (J7|) around q and keep 
terms up to 0(1) order. We obtain the equation derived 
by Escudero and Kamenev [H|]: 

y (Nw r + u r ) 



II. EIGENVALUE PROBLEM, WKB 
APPROXIMATION, AND FAST- AND 
SLOW-MODE SOLUTIONS 

When starting at t = from a sufficiently large pop- 
ulation, the probability distribution P n (t), as described 
by the master equation @ approaches, on a relaxation 
time scale t r , a long-lived metastable PDF peaked at a 
non-zero attracting fixed point of the rate equation. The 
metastable distribution is slowly "leaking" to zero, see 
Eq. (g]) . Let us denote the non-zero attracting fixed point 
by n = n* (n* = ri\^. for scenario A and B, respectively, 
see Figs. [Hand EH- Using Eq. (dJ), we arrive at an eigen- 
value problem for the QSD 7T n , n = 1, 2 . . . : 



[W r (n - r)7r„_ r - W r (n)TT n ] = -Eit n . (6) 



Importantly, the eigenvalue E = 1/r turns out to be 
exponentially small compared to the relaxation time t r . 
Therefore, the term in the right hand side of Eq. (|6]) can 
be neglected [HI, EH, HH, [22[ , and we have to deal with a 
quasi-stationary equation 



[W r (n - r)TT n ^ r - W r (n)-K n ] = , n = 1, 2, 



(7) 

For definiteness, we normalize the QSD to unity: 
Sn=i 7r ™ = 1- O nce 7r„ is found, we can use Eqs. (|3|) 
and g]) to calculate the MTE: 



E = 1/r = V W r (~r)n 



r<0 



(8) 



Let us introduce a rescaled coordinate q = n/N, where 
N 1 is the large parameter of the problem. The central 
assumption of our theory is that, after a proper rescaling 
of time which will be introduced shortly, the process rates 
can be represented as 

W r (n) = W r (Nq) = Nw r (q) + u r {q) + 0(1/N) , (9) 

where, for q = 0(1), w r (q) and u r (q) are 0(1). This 
assumption guarantees that the population be long-lived, 



-rS' 



1 



N 



S'i 



2N 



S" 



r w r 
N w, 



1 



0, (11) 



where the primes denote differentiation with respect to 
q. In the leading order O(N), this equation yields a sta- 
tionary Hamilton- Jacobi equation H (q, S') = 0, where 



F(?,p) = 5>r(ff) (e rp -l) 



(12) 



is the effective Hamiltonian, and p = S' is the momentum 
(12] ]. Therefore, in the leading WKB order, one needs to 
find zero-energy phase trajectories of the Hamiltonian 
(fT2j) . As w r (0) = for any r, one such trajectory is 
q = at an arbitrary p: the extinction line. This line 
is of no importance in the WKB theory, however. What 
we need are phase trajectories p = p{q). One of them 
is the relaxation trajectory p = 0. In general, there is 
one and only one additional phase trajectory for which 
p = p a (q) 7^ 0, except in some points q. Let us prove 
this statement. The Hamiltonian H(q,p) vanishes at p = 
0. Differentiating Eq. (Tl2^) twice with respect to p, we 
obtain H pp (q,p) = ^ Jr r 2 w r {q)e rp > 0. Therefore H is 
a convex function of p, and so it has one and only one 
additional real root [2J|. The relaxation trajectory p = 
and activation trajectory p — p a {q) give rise to the slow 
and fast WKB modes, respectively, as was shown in a 
particular example in Ref. [2~D ]. 

The q-dynamics along p — is described by the Hamil- 
ton's equation 



. dH(q,p) i 



H p (q,0) =^2rw r (q) 



(13) 



which is nothing but the (rescaled) rate equation (fTJ) . The 
nontrivial fixed points of the rate equation, qi = rii/N, 
are positive roots of the equation H p (qi,0) = 0. As 
Pa(qi) = 0, the activation trajectory p — p a (q) crosses the 
relaxation trajectory in these fixed points. As w r (0) = 
for any r, the small-g expansion of Eq. ()13|) generically 
starts with a linear term in q. In the remainder of this 
paper we assume that the linear decay rate is non-zero in 
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the leading order: w'_i(0) — a > 0. Therefore, we can al- 
ways rescale time, and all the rates, by a. This procedure 
uniquely defines the rescaling leading to Eq. ©. 

In extinction scenario A, see Fig. [21 the most probable 
path to extinction is the heteroclinic trajectory connect- 
ing the fixed point (g* , 0) of the Hamiltonian system (here 
g* is the attracting point of the rate equation) with the 
fluctuational extinction point (0, p/) . 

In extinction scenario B, see Fig. [31 the most probable 
path to extinction is composed of two segments. The first 
one is the activation trajectory: a non-zero-momentum 
heteroclinic trajectory connecting the fixed point (g*,0) 
with an intermediate fixed point (g rep ,0), where q rep is 
a repelling fixed point of the rate equation. The second 
one is a relaxation segment p = 0, connecting the point 
(g rep ,0) with the deterministic extinction point (0,0). 

For the fast mode one obtains UM 



S(q)=S^(q)= / d£p a (0, 



(14) 



where the integration constant is already accounted for 
by the prefactor A in Eq. (fTU]) . For the slow mode S(q) — 
S^{q) = 0. 

In the subleading 0(1) order, Eq. (fTTj) yields a first- 
order differential equation for Si(q): 



V w r e r P lrS[- r -p'- r -^)+ u r (e r P - 1) = , (15) 



where p = p a {q) for the fast mode, and p = for the slow 
mode. It is convenient to use the identities 

H q = J2 ™'MW P - 1) , ^ = E rw r ( q yp , 

r r 

H qq = <{lW P " 1) , H pp = r 2 Wr(q)e r P , 

r r 

H qp =Y,rw'Me rv , (16) 

r 

where the subscripts p and q stand for the partial deriva- 
tives. To remind the reader, all the rates in Eqs. (fTB"]) are 
rescaled with respect to the linear decay rate constant a. 
The fast-mode solution for S±(q) can be written as 

S[ f) (q) = J" dd 
H qp {LVa) + \H pp {i,p a )p' a {£,) - Z r M0(e rPa - 1) 



(17) 



where p a = p a (£) . This result was obtained by Escudero 
and Kamenev [221 ] in the context of stochastic population 
switches. The quantity H p (^,p a ) in the denominator of 
the integrand in Eq. (|17j) vanishes in every fixed point 
£ = q,; of the rate equation, including £ = 0. To see 
how the integrand behaves at the fixed points, consider 



the equation H[q,p a (q)] — for the activation trajectory. 
Differentiating it with respect to g, we obtain 

H q [q, Pa (q)]+H p [q,p a (q)}p' a (q) = 0. (18) 

One more differentiation gives 

H gq + H p p" a + (2H pq + H pp p' a )p' a = , (19) 

evaluated at p = p a (q)- By virtue of identities (fTB)) each 
of the first two terms of Eq. (fT9]) vanishes at q = qi and 
P = Pa{li) = 0, so the expression 2H pq + H pp p' a must 
also vanish there. As a result, the first two terms in 
the numerator of the integrand of Eq. (fTTj) cancel each 
other at £ = qi. The remaining term in the numerator, 
— 5Z r u r ({;)(e rPa ^ — 1), is proportional to (£— in the 
vicinity of £ = qi, exactly as the quantity H p {^,p a ) in the 
denominator. Therefore, the integrand is well behaved at 
£ = qi > 0. At q = 0, Si (<?) diverges logarithmically, see 
below. This divergence does not cause any concern, as 
the WKB approximation (which demands n > 1, or g > 
A^ 1 , for its validity) does not hold for small q anyway. 

One can partially perform the integration over q in 
Eq. ([T7) by using Eqs. (fT5|) and JT5J). After some algebra, 



where 



^ /) (g) = -ln v ^ 7 M + vi/( g ) : 



(20) 



2#„(£,Pa) 



Hp(^Pa) 



(21) 

and 5" (9) = p' a (q)- Note that 5" (9) does not vanish in 
any of the nontrivial fixed points g^, so the logarithmic 
term in Eq. ([20]) is well behaved there [H}. Indeed, from 
Eq. © S"(q) =p' a (q) = -H q [q,p a (q)}/H p [q,p a ( q )}. By 
Taylor-expanding the numerator and denominator in the 
vicinity of any nontrivial fixed point <?,•, one can see that 
S" (Qi) 7^0. 

Therefore, the general fast-mode solution is given by 
Eq. (0 with S (f Hq) from Eq. (UJ), and ^(g) from 
Eqs. (f20]) and ([21]) . It also includes a constant prefac- 
tor A — Af which can be found immediately. Indeed, 
at W > 1 the QSD is strongly peaked around the at- 
tracting fixed point q — q*. Here the fast-mode solution 
dominates, and Af can be found by normalizing to unity 
the gaussian asymptote of the QSD around q = q*. The 
gaussian asymptote is obtained by expanding the QSD 
pop in the vicinity of q = g*: 



n(q) ~ A f t 



-NS if) (q«) 



) (q,)-(N/2)S"( q ,)(q- q ,) 2 



(22) 



where we have used the equalities S'(q*) = p a (<l*) = 0. 
To see that S"(q*) > 0, one can again use Eq. (TKfl) . At 
q = g* it reads 2i/ p9 (g»,0) + H pp (q t ,,0)p' a (q t ,) = 0. By 
virtue of Eq. (fT6)) H pp (q*,0) > 0. The quantity H pq (q*,0) 
is the g-derivative, evaluated at g = g», of the expression 
in the right hand side of the rate equation (flU]) . As the 
point q — g* is by assumption attracting, H pq (q*, 0) < 0. 



() 



Therefore, p' a {q*) = S"(q*) > 0, and the asymptote (|2"2"|) 
is indeed a gaussian distribution. Normalizing it to unity, 
we obtain 



A, 



S"{<1*) A r 5 (/) (g»)+5 1 (/) (g») 



2nN 



so the fast-mode solution is fully determined: 



w(q) 



S"{q*) ms<n (q*)-s<-n (g)]+s(« ( q ,)-s[ f> ( q ) 
2irN 



with SW(q) from Eq. (HJ) and S^fa) from Eqs. O 



(23) 



,(24) 



and (|2Tl) . 

Now consider the slow-mode solution for which 
S( s \q) = 0. The subleading-order contribution s[ s \q) 
is found by putting p = in Eq. (fT7|) : 



Si'Hq) 



hxH p (q,0). 



(25) 



Then Eq. (1101) yields the general slow-mode solution for 
the QSD 21~E|: 



^s(q) 



H p {q,QY 



(26) 



where A s is an arbitrary constant. The minus sign is put 
here for convenience, because in the region of < q < q\, 
where the slow- mode solution is relevant (see Section PVT) . 
q = H p (q, 0) < and A s > . One can see from Eq. (j26|) 
that the slow-mode solution diverges in the fixed points 
of the rate equation [2l[. This divergence will be cured 
in Section [V] 

We show in the following that, for a given extinction 
scenario and in a given region of q, only one of the modes, 
either fast or slow, dominates the resulting QSD, while 
the other one must be discarded. Before we deal with this 
issue, however, we recall that the WKB approximation 
breaks down at n = 0(1). To find the QSD for all n we 
will solve Eq. (J7J) in the region of n <C N by recursion 
and then match the recursive solution with either the 
fast- mode (in scenario A), or the slow- mode (in scenario 
B) WKB solution in the joint region of their validity. 



where only processes with w' r (0) ^ contribute. One can 
look for particular solutions of this recursive equation in 
the form 7r„ = f n /n thus arriving at an equation with 
n- independent coefficients: 



X>r(0)(/n-r-/n)=0. 



(28) 



In the remainder of this paper we make the following 
simplifying assumption: 



,(0) = <,(0) = •••=(). 



(29) 



That is, we assume that the rates of the multi-step loss 
processes n — > n — m, where m = 2, 3, . . . , do not have, 
in the leading order in JV, linear terms in their Taylor 
expansion in n. This assumption is always satisfied for 
stochastic chemical reactions (where pairs, triplets, . . . , 
of reacting particles are needed to bring down the number 
of particles by 2, 3, ... ). The conditions ([29jl also hold for 
all models of population biology and epidemiology we are 
aware of. 

Using Eq. (|29|) and the equality w^ 1 (0) = 1 (to remind 
the reader, we are using rescaled variables), we rewrite 
the recursive Eq. (|2"5|) as 



fn+l — 



1 



K 

£ 

r=l 



w' r (0) 



fit 



K 

£* 

r=l 



(30) 



where K = r max <C N. If there is no degeneracy, the 
general solution of Eq. (|3"0|) is a linear combination of 
all particular solutions /„ = A~™, where A obeys the 
characteristic polynomial equation of degree K + 1: 



K 



w ;(o)A r+1 - 



i 



£>;(o) 



A + 1 = 0. (31) 



Note, that A = Ao = 1 is always a root. Let us show that 
Eq. (I31[) has one and only one additional positive root, Ai, 
while all others roots A2, A3, . . . , Ak are either negative 
or complex. First, we establish a connection between 
the roots Ai of Eq. pip and the crossing points with 
the p-axis of the zero-energy trajectories of the WKB 
Hamiltonian (JTSJ) . By expanding w r (q —> 0) ~ w' r (Q)q, 
Eq. ([12"]) becomes 



X(0) (e"-l)=0, 



(32) 



III. RECURSIVE SOLUTION 

The objective of this Section is to approximately solve 
Eq. ([7|) at sufficiently small n. The exact criterion of 
smallness will appear later, when we match different so- 
lutions in joint regions of their validity. 

In the leading order in N we take W r (n) = Nw r (n/N), 
see Eq. ©, and expand it in n/N up to the linear term: 
W r {n) ~ Nw r (0) + nw' r (0) = nw' r (0). Then Eq. 
becomes 



4.(0) [(n - r)7r„_ r - nir n ] = , 



(27) 



where only terms with w' r (0) 7^ contribute. Putting 
e p = A and ui'_ 1 (0) = 1 and using Eq. (f2T))) . one can see 
that Eq. (I3"2"j) coincides with Eq. (f3"Tj) . As we have shown 
that the equation H (q, p) — has, for any q, two and only 
two real solutions for p [23J], Eq. (j3T|) also has two and 
only two real solutions for A, both of them positive. The 
roots Ao = 1 and Ai correspond to the crossing points 
with the p-axis of the slow and fast modes, respectively. 

Dividing Eq. (|3ip by 1 — A, we arrive at a polynomial 
equation of degree K: 



K 



£<(0)(A 



A' 



0. 



(33) 
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For K < 4 the roots of this polynomial can be expressed 
in radicals. For K > 5, they need to be computed nu- 
merically. Assume that we have found all of the roots of 
Eq. ([55)1 Aj, i — 1, 2, . . . , K. If there is no degeneracy, 
the general solution for /„ is 



K 



K 



fn — " — Co + ^ C l \ i n , 



(34) 



i=0 



i=l 



where the coefficients Ci are the following (see Appendix 
A for the derivation): 



(35) 



n ( a »-^) 

3 + i 



The coefficient Co, corresponding to the root Ao = 1, 
can be expressed through the coefficients w' r (0), r = 
1,2, ... ,K (see Appendix A): 



C = 



h 



1 - tui(0) -2w' 2 (Q) Kw' K (0) 



(36) 



Before writing down the general solution of the recur- 
sive equation (|27p for the QSD, we recall a simple relation 
[25l ] between /1 and the MTE t. Using the rescaled re- 
action rates, we can rewrite Eq. (jSJ as 

T^ 1 = a 2_. W r ( — r) 7T_ r . 
r<Q 

In view of the conditions (|29p . only one term in the sum 
survives in the leading order. As v/_i(0) = 1, we obtain 

t -1 = aW_i(l)7ri ~ aiTi — afi , 

so fi ~ l/(ar). Now we switch to the rescaled variable 
q = n/N, use the relation ir(q) = f n /(Nq) and Eq. AMI) , 
and obtain the small-q asymptote of n(q): 



K 



n(q) 



j'=i 



E 



A 



— JVqr 



(37) 



1=0 n (^-a.) 
3=0 

3 + * 

The validity region of this asymptote (which includes the 
yet unknown t to be found later) is scenario-dependent. 
It is a relatively narrow region q <C N^ 1 / 2 in scenario 
A, and a broader region q <C 1 in scenario B. The differ- 
ence comes from the fact that in scenario A the recursive 
solution needs to be matched, at n 3> 1, with a rapidly 
growing fast-mode solution, whereas in scenario B the 



matching needs to be done with a slowly varying slow- 
mode solution, see Sections |IV] and |Vl respectively. 

What is the role of complex roots of the polynomial 
equation ((33]) in the recursive solution ((37)) ? These can 
appear only for K > 3, and they come in complex con- 
jugate pairs: Xj and A& = Aj. One can show, by using 
Eq. (I3"5j) . that the coefficients Cj and Ck, corresponding 
to Aj and Aj, are also complex conjugate: Ck — Cj, so 
that n(q) from Eq. ((3"T|) is real-valued as expected. When 
complex roots are present, the QSD at small n may ex- 
hibit rapidly decaying oscillations as a function of n. 

Let us now determine then3>l, org^>A _1 , asymp- 
tote of the QSD ((3"T|) . in each of the two extinction sce- 
narios. This asymptote will be matched, in each scenario, 
with the dominant WKB mode. Equation (|37p includes 
K + 1 terms. At n 3> 1 the leading contribution comes 
from the term with the smallest |Aj|. The rest of the 
terms are exponentially small compared to the leading 
one and can be safely neglected. 

In scenario A, the two positive roots of Eq. (f3"Tj) are 
< Ai < 1 and Ao = 1, whereas the rest of the (negative 
or complex) roots obey the inequality |A<>i| > Ai, see 
Appendix B. In this case the asymptote of the recursive 
solution (|3"T)) at n ^S> 1, or q 3> TV -1 , is 



Tr(q) 



Nq 



(38) 



arNq ' 

where the positive constant A\ (see Appendix B) satisfies 



A 1 



(39) 



(-l)*fU- 

j=i 

K 

n ( a i-^) 
.7=0 



In this case Ai = e Pf corresponds, in the WKB-language, 
to the pf < crossing point of the activation trajectory 
and the p-axis, see Fig. [2] Therefore, to set the ground 
for matching Eq. (f3"81) with the WKB solution in scenario 
A, we can rewrite Eq. (f3"5|) as 



A 1 e -Npsq 
arN q 



(40) 



In scenario B the n>l asymptote of Eq. (|57jl is quite 
different. Here the root of Eq. (|3"T|) with the smallest 
absolute value is Ao = 1, see Appendix C. Therefore, the 
i = term in Eq. (|37|) is dominant, and we obtain 



n(q) 



1 



arNq[l - w[(0) 



2*4(0) Kw f K (0)] ' 

(41) 

where we have used Eqs. (|35|) and ([36]) . The asymp- 
tote (|41|) can be expressed in terms of the WKB Hamil- 
tonian ([T2]). Using Eq. JIBJ, we obtain H qp (0,0) = 



J2r rw r(Q)- Recalling that w'_ 1 (Q) = 1 and using 
Eq. J29]), we can rewrite H qp (0,0) as 

H qp (0, 0) = Kw' K (0) + (K - 1)^(0) + • • • + w[ (0) - 1 . 

(42) 

As q = is an attracting point here, H qp (0, 0) < 0. Then, 
using Eq. (j42|) . the asymptote (|4T|) becomes 



where 



?r(g) 



tf|# w (0,0)|g 



(43) 



Note that Ao = 1 = e Ps corresponds to the zero- 
momentum (p s = 0) crossing point of the relaxation tra- 
jectory and the p-axis, see Fig. [3J 



IV. EXTINCTION SCENARIO A 
A. General case: multi-step processes 

In this section we calculate the MTE and QSD for ex- 
tinction scenario A. Here extinction occurs along the acti- 
vation trajectory: the heteroclinic trajectory, connecting 
the metastable point (n\, 0) and the fluctuational extinc- 
tion point (0,Pf) of the phase plane (n,p), see Fig. [2] 
In this case the slow-mode solution is negligible com- 
pared to the fast-mode solution in the entire region of 
q > 0. Furthermore, the fast-mode solution ([2"4"|) can 
be directly matched with the recursive solution (|37j) in 
the joint region of their validity which turns out to be 
1 < n < A 1 / 2 , or N^ 1 < q < A" 1 / 2 . 

To implement the matching procedure, we first find 
the q <C TV -1 / 2 asymptote of the fast-mode solution ([24)) . 

Because of the divergence of S\ J '(q) at g = 0, we should 
proceed with care. Let us rewrite Eq. as 



, , = \JS"{qi) p N[S i -f\ qi )~S^\q)\ + S[ f) (qi)-[S[ S) (q)-\nq\ 

V2^Nq 

Here we have introduced the 1 /q prefactor which diverges 
at q = 0, and made up for it by adding lng in the expo- 
ncnt. Let us show that the expression S\ (q) — \nq in 
the exponent is regular at q = 0. We represent hiq as 
J 9 d^/C and use Eq. (JTTJ) to rewrite S[^\q) — Inq as an 
integral over £. Now we Taylor-expand the integrand 

H pq (t, Pa ) + \H pp {tp a )p' a {0 - £,«r(£)(e rp " - 1) 1 



Hp(i,Pa) 



in the vicinity of £ = up to linear terms. The divergent 
terms cancel out, and the remaining expression 



(l/2)H ppq (0,p f )p' a (0) - EX(0)(e^ - 1) 

H pq (0,pf) 

is finite. Now we rewrite Eq. (|44l as 



. (45) 



, , = \JS"{qi) gi p JV[s»)( gl )-S(»(g)]+0( gl )-0(g) / 46 n 
^/2ttN q 



4>(q) = S\ f \q)-\nq 



(47) 



is regular at q = 0. By Taylor-expanding the exponent 
of Eq. (|46|) around q = to first order, we obtain the 
q <C A -1 / 2 asymptote of the fast-mode solution: 



V2ttN q 



(48) 



This asymptote can be matched with the asymptote of 
the recursive solution at ./V -1 <C q <C A r_1 / 2 , given by 
Eq. ([40]) . This matching yields 



A 1 V2tt 



,JV[S(«{0)-SW( 91 )]+«0)-*(ft) 



aq iy /NS"( qi ) 



(49) 



where </>(g) is given by Eq. (I47|) . a is the linear decay rate 
constant in physical units, and A x is given by Eq. ((3"9")) 
(26| . The general expression ([4^]) for the MTE in scenario 
A is one of the main results of this work. The leading 
term in the exponent, proportional to A, is the effective 
entropy barrier to extinction. The proportionality factor 
is the absolute value of the area under the activation 
trajectory, see an example in Fig. [5] 13]. Noticeable is 
the presence of the large factor A 1 / 2 in the pre-exponent. 
The constant has a clearly non-WKB nature, as it 
comes from the recursive solution of the quasi-stationary 
master equation at small n and is contributed to by all 
of the roots Xi, i = 0, 1, . . . , K. 

Another important result is the QSD in extinction sce- 
nario A. It is determined by the asymptotes (|2~4")) and ((3"?]) 
which coincide, in the leading order, in their joint region 
of validity A -1 < q < A~ 1/2 , or 1 < n < A 1/2 . 



B. Single-step processes 

Remaining within scenario A, we now turn to an 
important sub-class of stochastic population processes: 
single-step processes. Here there are only two non-zero 
process rates: W±i(Nq) = W±(Nq) = Nw±(q)+u±(q) + 
. . . , where all the rates are normalized by the linear de- 
cay rate constant w'_(0). In this case the expressions for 
the MTE and QSD can be simplified considerably. The 
WKB Hamiltonian (TT21 becomes 

H(q,p) = w + (q)(eP - 1) + w_(q)(e-r - 1) . (50) 

The rate equation is q = w + (q) — u>-(q). In scenario A 
one has w' + (0) > w'_(0) = 1. Here it is convenient to 
denote the ratio of the linear birth and death rates by 
R ee w' + (0)/w'_(0) = w' + (0). For R > 1 the fixed point 
q = of the rate equation is repelling. The activation 
trajectory is p a {q) = — ^ D -[w+(q)/w-(q)}, and 



S U \q) 



In -— d£. 



(51) 
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Now we calculate the following quantities on the activa- 
tion trajectory: 



Pa(0) = S "(l) = ± > HpiliPa) =W--W + , 

W-w' + u>+w'_ 
Hp q {q,Pa) = , H pp (q,p a )=w-+w + , 



--u, I — -!)+«_ | — - 1 I .(52) 



Substituting these into Eq. (fTT|) , we obtain after simpli- 
fications 



r 




_ u + \ 






W+J 



<% + ilnK (q)w- (<?)]. (53) 



Plugging this into Eq. (|4"?| yields 



w+(q)w-(q) 



(54) 



Now we can calculate e^ )-"^) which enters Eqs. (gSJ) 
and ggj) : 



o 0(O)-0(gi) 



<7i 



w+(<7i 











exp 








Jo \ w + 







(55) 



where we have used the following relations: (i) as q — > 
0, W+ (q)w^(q)/q 2 -> <(())«/ (0) = fl, (ii) w+fa) = 
w -(qi) ( as 9i is a fixed point of the rate equation), and 
(iii) w'_(Q) — 1 because of the rescaling of the rates. 

Now we turn to the recursive solution at small n, pre- 
sented in Section HTT1 For single-step processes K = 1, 
and so Eq. (f3"H has only two roots: X = 1 and Ai = 
l/if+(0) = 1/R. Therefore, rewriting the small-g asymp- 
tote ([37]) of the QSD in terms of 7r„ , we obtain 



(R-l)n 



whereas the constant A\ from Eq. ((39, 
ging this constant and Eqs. (f5Tj) and 
the MTE, we obtain 



V2nR 



)dg 



a(R-l)w+(q 1 )y/NS"(q 1 ) 



(56) 

is 1/(R-1). Plug- 
SS) in Eq. (ggj) for 



(57) 



In the particular case u± = u_ = Eq. (|5T|) coincides 
with Eq. (19) of Ref. [la], obtained, via a saddle-point 
approximation, from the exact expression for the MTE 
of a single-step process. Doering et al. assumed in their 
derivation that the subleading contributions u± to the 
process rates W±i(ri), see Eq. (j9|), vanish. As a result, 
the factor exp [J^ 1 (u + /w + — w_/w_) dq\ is absent from 
their Eq. (19). While the assumption u + = «_ = may 
hold in some simple models, it does not hold in general. 
For example, it does not hold for stochastic chemical re- 
actions where the rates are combinatorial, as in one of 
the examples we present in subsection IIVD below. 



C. Extinction near transcritical bifurcation point 

Now let us return to a general set of (not necessarily 
single-step) processes. Our objective is to simplify the 
MTE (|4"9"|) in the special regime when the population, 
as described by the rate equation, is very close to the 
characteristic (transcritical) bifurcation point of scenario 
A. Here the attracting point q = qi is very close to the 
repelling point q = 0, so that q\ <C 1. This also implies 
\pf\ <§; 1 [3, [27fl. Taylor-expanding Eq. (fT2"|) in q and p 
around q = p = 0, we obtain 



H(q,p)^qp^2 [«4(0) 



= 0. 



|r<(0) + |rV(0) 

r 

The trivial solutions are the extinction line q = and 
the relaxation trajectory p = 0, whereas the nontriv- 
ial solution yields a straight-line activation trajectory. 
Using Eq. (fTrj|) and expanding the algebraic equations 
H p (q,0) = for qi, and H q (0,p) = for pf at small q 
and p, we can represent the activation trajectory as 



Pa(q) 



1 



Here 



Qi 



2g gP (0,0) 
H qqp (0,0) ' 



where H qqp (0,0) < 0, and 



Pf 



2H qp (0,0) 
H qpp (0,0) ' 



(59) 



(60) 



(61) 



where H qpp (0, 0) 



> 0. Exactly at the bifurcation the 
rate constants are such that H qp (0,0) = 0. Here the 
attracting fixed point q\ merges with the repelling point 
q = 0. The coordinate of the attracting fixed point qi = S 
can serve here as the distance to the bifurcation. [The 
third derivatives of the Hamiltonian, which appear in the 
denominators of Eqs. f|60|) and (|6"T|). are generically of 
order unity] 

Now, using Eq. (|5"9")1 , we can calculate S"(qi) = 
Paili) = ~Pf/qi — ~Pf/$ > 0; the fast- mode action 



S (J) (q) 



Pa{0d£,=Pfq 



8 2 



(62) 



and the accumulated action between the points q = q\ = 
S and q = 



AS = S {f) (0)~S {f) (q 1 ) = - 



Pfli 



PfS_ 
2 



> . (63) 



This quantity is the area of a triangle [l4( , see Fig. 2J 

To find the fast-mode correction to the action, S{ (q) , 
we expand Eq. ([T5|) in the vicinity of q = and p = 0, 
keeping only the leading order terms. This yields 



(64) 
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ft, 





q 



FIG. 4: (color online). Shown are typical zero-energy tra- 
jectories of the WKB Hamiltonian (|12[l in scenario A, close 
to the bifurcation where q\ — 8 <C 1, and pf <C 1. Here, 
the activation trajectory is a straight line, and AS given by 
Eq. (|63[) is the area of the shaded triangle. 



Although Eq. (JTHJ) breaks down at 6 ~ N^ 1 / 2 , one can 
still predict a scaling relation for the MTE in this region: 
r ~ N 1 / 2 /a. The symbol ~, here and in the following, 
means "of the same order as" . 



D. Examples 

We will now illustrate our theory by calculating the 
MTE in four pedagogical examples of extinction scenario 
A. The first three of them are single-step processes: the 
logistic Verhulst model of population dynamics, a set of 
three chemical reactions, and the SIS model of epidemics. 
The fourth example - another set of three chemical reac- 
tions - involves a two-step process. We will also consider 
all of these examples near the bifurcation. 



whereas the subleading terms u r in the rate expansion 
do not contribute. The solution of this equation is 

S[ f \q) = \nq. Then, by virtue of Eq. g7|), we obtain 
4>{q) = 0. That is, near the bifuraction, the subleading 
WKB correction vanishes. 

Now let us consider the coefficient A x [see Eq. ([M]) ] 
which enters Eq. ([4"9"|) . Among the roots \ of the polyno- 
mial (|31[) , which contribute to A\ , two are special: Ao = 1 
and Ai = e PS . Near the bifurcation \pf\ <C 1, so we can 
write Ai ~ 1 + pf < 1. As a result, A\ = Ai/pf. where 



c-i)*iU 

3=2 

J=2 



(65) 



is a negative constant of order unity, and we have put 
Ai = 1 in the expression for A±. Furthermore, the roots 
Aq, A3, ... , Xr of the polynomial (j33|) can be evaluated 
at the bifurcation point. As a result, one can express Ai 
via the linear branching rates ui' r (0). After some algebra 



Ai 



K 



5>(r + lK(0) 



(66) 



Substituting all of the above into Eq. (|4"9"|) . we obtain the 
MTE close to the bifurcation point: 



2tt \A x \ 



N aDl S 2 6Xp V 



(ND 2 5 2 



(67) 



where and D A = ^J\H q(jp (0, 0)\/H qpp (0, 0) = 0(1) 
should be evaluated at the bifurcation. Equation l|67p is 
valid when A^ 2 ^> 1. For sufficiently large N this strong 
inequality is compatible with the strong inequality 8 <C 1 
which describes closeness to the bifurcation. Note that 
the constant A\ = 0(1) is determined by the full small-n 
recursive solution that we found in section IIIII 



1. Verhulst model 

The generalized Verhulst model is a stochastic logistic 
model: a single-step Markov process with birth and death 
rates 



W+i = W+ = ot\n — a-in 
W-! = W- = p in + p 2 n 2 , 



(68) 



respectively, where ai,a.2,0i and j3 2 are non-negative 
rate constants. The quadratic corrections account for 
competition for resources [28[. It is customary to put 
Q2 = in Eq. (|68|) fl8j ] . and this is what we will do here. 
Rescaling time by the linear death rate constant we 
bring the rates to the form given by Eq. ©: W+ = Bn 
and W- = n + Bn 2 /N, where B — cti/fti is the ratio 
of the linear birth and death rates, and N = B/3i/^2- 
According to Eq. © w + = Bq, W- = q + Bq 2 , and 
u + = U- = 0. At B > 1 the fixed point q = of the rate 
equation is repelling, whereas qi = 1 — l/_B>0is attract- 
ing. Here we have S"{q{) = p' a {q x ) = 1, A x = l/(B - 1) 
and 



In dq 



1-B + lnB 



/o w + B 
Therefore, the MTE [Eq. (|57p] in physical time units is 



1 [2^ 4b 



fa V N (B - I) 2 



exp 



N 



B-1-biB 
B 



, (69) 



which coincides with previous results obtained by differ- 
ent methods [H, HI] . 

In this simple example the process rates satisfied the 
conditions u + = w_ = 0, so Eq. (|6"9"]) could have been 
obtained from Eq. (19) of Ref. [18|]. In the next example 
we relax one of the two conditions and show that, as 
predicted by our more general Eq. (|57[) . the pre-exponent 
of the MTE changes. 
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2. Reactions A 2A and A 







Consider a set of three reactions among particles A: 
branching A A 2A, a reverse reaction 2A A A, and de- 
cay A A 0. As observed in Ref. [27|, this set of reactions 
can be viewed as a generalization of the Verhulst model 
considered in the previous example. Indeed, by impos- 
ing a special relation, a — 2(/x — 1), between the rate 
constants, and by denoting fi = 1 + B/N and A = B, 
one recovers the process rates w± and it± of the Verhulst 
model. For this special choice of rate constants one has 
ii + = w_ = which yields Eq. ([Ml for the MTE. Now, 
what if the rate constants a and /i are independent? As 
usual, we can normalize time and the reaction rates by 
fi and denote B = X/fi and N = 2X/a. By virtue of 



Eq. © we can write w + — Bq, w_ = q + Bq 



' + 



= 

and U- = Bq. The rescaled rate constants are identical 
to those of the Verhulst model, except that U- is now 
nonzero. As a result, 



exp 



and Eq. ([57]) for the MTE yields 




exp 



a \ 



dq 



B. 



A 



(70) 



where we have returned to physical units. This result 
cannot be obtained from Eq. (19) of Ref. [18j. 



Therefore, the MTE [Eq. (jSTj) ]. in physical time units, is 
given by 



1 2n Rn 



n v n (R - iy 



exp 



-V ( lni? + -^-l 
H 



,(71) 



which coincides with previous results obtained by differ- 
ent methods (H-tH. 



4- Branching- annihilation- decay 



Now we consider another set of stochastic reactions 
among particles A which include, in addition to single- 
step processes A A 2A and A A 0, a two-step process: 
binary annihilation 2 A A- 0. This pro blem was previ- 
ously solved by Kessler and Shnerb [16j. Here we show 
that their result for the MTE follows from our Eq. 09} . 

In our notation, the transition rates between the states 
n and n + r are given by 

anln — 1) . , 

Wi = Xn, W-i = [in , and W^ 2 = — -z ■ (72) 

Rescaling time fit — > t and denoting Ro = X/ /u, and N = 
X/a, we obtain Eq. ([9|) with 



«h = R q , tw_i = q , = 



R Q q 2 



u\ = , w_i = , u_2 = — 



R q 



(73) 



3. SIS model 

Now let us consider the well-known SIS model of epi- 
demics, see e.g. Refs. (2^ . |29| and references therein. The 
SIS model deals with dynamics of a population which 
consists of two groups of individuals: susceptible to in- 
fection and infected. It is assumed that infection does 
not confer any long-lasting immunity, and infected indi- 
viduals become susceptible again after infection. When 
demography (births and deaths) is negligible, the to- 
tal number N of individuals in the two groups is con- 
served. As a result, the model becomes effectively single- 
population, with the effective rates 



W+ = Xn{N - n) , W. 



fin . 



Mathematically, this model is just another example of 
the generalized Verhulst model, see Eq. (|68|) . where one 
chooses /3i ^ but /3 2 = HI- 

Let us denote Ro = XN//I and rescale time and rates 
by the linear decay rate constant /i = XN/Rq. The 
rescaled rates become u>+ = Ro(q — q 2 ) and w_ = q, 
while u + = u_ = 0. The fixed point q\ = 1 — l/i?o of the 
rate equation is attracting when Ro > 1. Furthermore, 
S"(qi) = p' a (qi) = Ro, and A x = 1/(Rq ~ 1) (as in the 
above notation R = Rq). Finally, 



In dq 

w + 



1 



1 



lni? . 



In the rescaled notation, the attracting fixed point is q\ = 
1 — 1/Rq which demands Ro > 1. The WKB Hamiltonian 
(ITS]) takes the form 

H(q,p) = R q (e* - 1) + q (e^ - l) + ^ (e" 2 * - l) . 

(74) 

Solving the equation H[q,p a (q)] = 0, we obtain the acti- 
vation trajectory 



(75) 



where u = 2 + qRo and v = \Ju 2 + 8gi?Q. The zero- 
energy phase trajectories of this system are shown in 
Fig. El 

Now we use Eqs. JUJ), ([T7]), and (@7J| and obtain 



AS = S^ n (0)~S u ' ) (q 1 ) 



2 1- 



and 



Ro 



In 



2 V + Ro 



2R 



^(Ro + l)(3R -l) 



(76) 



(77) 



Furthermore, S"(qi) = 2R /(3Ro — 1) and, as in the 
previous examples, the only root of Eq. ([55)1 is Ai = 1/Rq- 
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Therefore, by using Eq. (13T)[) . we have A\ — 1/(Rq — 1). 
Substituting all of the above into Eq. (|49|) . we obtain, in 
physical time units, 



2^ 



R, 



3/2 



which coincides with the result of Ref. UM ■ 



5. Examples 1-4 near the bifurcation point 

Because of their simplicity the examples 1-4, presented 
above, give identical results for the MTE near their cor- 
responding bifurcation points, described by the equation 
H qp (0,0) = 0. The small distance to the bifurcation 5 
in all these examples is the ratio of the linear birth and 
death rates minus 1: B — 1 in example 1, X/fi — 1 in 
example 2, and Rq — 1 in examples 3 and 4. In all four 
examples Da = y/\H qqp (0, 0)\/H qpp (0 : 0) = 1. As there 
is only one linear branching process in each example, one 
has K — 1, so Eq. ([5FJ|) yields A\ = — 1. As a result, in 
all four examples 



/2tt 



[N5< 



■ ex P i 
NS 2 V 2 



(79) 



where a denotes, in each example, the linear decay rate 
constant. To remind the reader, Eq. ([75)) is valid when 
N5 2 S> 1 which, together with 6 <C 1, yields the double 
inequality N^ 1 / 2 < 6 < 1. At 5 ~ N^ 1 / 2 Eq. ([79]) 
predicts the following scaling relation for the MTE: r ~ 
N 1 / 2 /a, in agreement with Ref. [281 ]. 



V. EXTINCTION SCENARIO B 
A. General case: multi-step processes 

In this section we calculate the MTE and QSD for ex- 
tinction scenario B. Here extinction occurs along a trajec- 
tory composed of two segments: the non-zero-momentum 
heteroclinic trajectory connecting the hyperbolic fixed 
points (^2,0) and (nj.,0) (the activation trajectory), and 
the zero-momentum segment going from n = ri\ to n = 
(the relaxation trajectory), see Fig. [3l 

A straightforward way to calculate the MTE starts 
with finding the WKB solution for the QSD at n > 1. 
Then one should match it with the small-n recursive so- 
lution (I37p . as in scenario A. The matching region in this 
case is 1 <C n <C N, or N^ 1 -C q <C 1. After having 
found the QSD, one can determine the MTE by using 
Eq. (O. 

Actually, there is a shortcut to finding the MTE which 
does not require the knowledge of the small-n recursive 
solution. This is because the solution includes a constant 
probability current flowing from a close vicinity of n = 
m to a close vicinity of n — 0, as shown below. This 



probability current is equal to the escape rate from the 
metastable state q = q2, and it is determined by the 
WKB-asymptote of the QSD, with no use of the small- 
n recursive solution. One of the objectives of our work, 
however, is to also find the QSD of the metastable state, 
and therefore we will follow the straightforward way. 

In contrast to scenario A, where the WKB solution is 
determined solely by the fast mode, in scenario B the 
WKB solution is more complicated. Here the fast mode 
dominates to the right of the point q = qi (but not too 
close to qi), whereas the slow- mode solution ([2"fi[) dom- 
inates at < q < qi (again, not too close to q = qi). 
Furthermore, the slow- mode solution diverges at q = qi, 
and curing this divergence demands going beyond the 
WKB approximation in a boundary layer \q — q±\ <C 1 
where the fast and slow modes are strongly coupled. As 
a result, the QSD at n > 1 involves three distinct asymp- 
totes which need to be matched to one another. All this 
is very similar to what happens in other types of popu- 
lation escape problems: to an absorbing state at infinite 
population size [2l| or to another metastable state [HI]. 
Much of the calculation is very similar to that of Refs. 
[HI [HI; but we will present it here for completeness. 

In the boundary layer \q— qi\ <C 1 the momentum p is 
small, that is fluctuations are weak. Here we can apply 
the van Kampen system size expansion 0] to the quasi- 
stationary master equation ([7]). Let us denote f(q) = 
W r (q)n(q) ~ Nw r (q)Tr(q) [it suffices to keep only the 
leading term in Eq. (J5J)]. Taylor-expanding f(q — r/N) 
around r = 0, we obtain 



f(q r/N) ~ f(q) - ^f(q) + ^f'(q) ■ (80) 



Plugging Eq. ([80[) into Eq. ([7|) and integrating once, we 
obtain 



where J = const. Now, 



(81) 



f'(q)=N{7T'(q)w r (q)+Tr(q)w' r (q)}, 

but n'(q) ~ Nn(q), so the second term in f'(q) is negli- 
gible. Therefore, we obtain 

2 

- 7r(g) rw r(<l) + A?) E ^N Wr{q) = J - (82) 

r r 

The first term on the left corresponds to drift, the sec- 
ond one to diffusion. With the diffusion neglected, one 
obtains a (slow-mode) solution for ir(q) which diverges 
at fixed points of the rate equation. The diffusion term 
cures this divergence by providing coupling between the 
slow mode and fast modes, as observed in Ref. [21[. 

Now we use Eq. (fTH)) and evaluate the drift and diffu- 
sion terms in the vicinity of q = q%. In the drift term 



^2rw r (q) = H p (q,0) 



{q- q 1 )H pq (q 1 ,Q) 
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while in the diffusion term it suffices to put w r (q) ~ 
w r (qi). Denoting x — (q — qx)/l, and I 2 = 
H pp (qi,0)/[NH pq (qi,0)] [as one can check, both 
H pp (qi,0) and H pq (qi,0) are positive], we obtain the 
boundary-layer equation [2l| 



tt'(x) — 2xn(x) = J , 



(83) 



where the rescaled constant current J is to be found later. 
The general solution of Eq. ([53"]) is 



ir(x) = c\e x2 + ^erf(x) . 



(84) 



where c\ is another constant. Now we can match this 
solution to the slow- mode solution at x < and \x\ 3> 1, 
that is, at N^ 1 / 2 -C qi — q -C 1. To eliminate the expo- 
nential growth at x < 0, one must choose ci = J^/tt/2, 
so the asymptote of the boundary-layer solution l|84p at 
— x 3> 1 becomes 



7r(a;) 



J 



Hp P (qi,0) 



J 



2x 2( 9 i-g)V NH pp ( qi ,0) 



(85) 



The slow-mode solution (|26|) at g ~ gi can be approxi- 
mated as 



-4., 



Matching the two asymptotes, one obtains 

2A,VN 

J 



^/H p p( qi ,0)Hp q (q u 0) 



(86) 



(87) 



To find the still unknown constant A s , we have to match 
the x ^> 1 asymptote of the boundary-layer solution 
which is 



7r(x) 



^Hpp( qi ,0)H pq ( qi ,0) 



x exp 



NH pq ( qi ,0) 
Hpp(qi,0) 



(q-qi) 2 



with the asymptote of the fast-mode solution at N 1 / 2 <C 
q — qi <C 1, which is 



n(q) 



x e 



S"(q 2 ) 
2ttN 

N[S^(q2)-S ( fHq 1 )]+S[ f> (q 2 )~S l 1 f) (q 1 )-{N/2)S"(q 1 )( q - qi ) 2 



(89) 



Here we have used the equalities S'(qi) — p a (qi) = 
and neglected terms of order q — qi <C 1 in the exponent. 
Putting q = qi into Eq. (fit?)), we obtain 



p' a (qi) = s"( qi ) 



2g M (gi,0) 
H PP {qi,0) ' 



(90) 



where S"(qi) < 0. Matching the asymptotes 
([89]) and using Eq. ([20]), we find 



and 



A s = 



Hpp(q ly 0)^/\S"(q 1 )\S"(q 2 ) 
AttN 

„N[S<«( 92 )-S(«( ?1 )] + [S< /) (<, 2 )-S< /) (g 1 )] 



(91) 



What is left is to find the MTE by matching the slow- 
mode solution at q <§; 1 with the recursive solution (|4"3")l 
at q > N^ 1 . Using Eq. ([25]) we obtain, at q < 1: 



7r s (g) 



-4., 



gi? P9 (0,0) g|ff P9 (0,0)| 



(92) 



where H pq (0, 0) < is given by Eq. ([42]) . Comparing this 
with Eq. (|4*3^) and using Eq. ([{?!"]) ■ we obtain 



47T 



^ pp (g 1 ,0) x /|^'(g 1 )|^"(g 2 ) 



x e JV[s(«( 9l )-s(«(« 3 )]+[sF ) (9i)-sF ) (?2)] ) (93) 

where a in the linear decay rate constant in physical units 
[26j . The expression ([M]) for the MTE in scenario B is 
an important result of our work. The leading term in 
the exponent, proportional to N, is the effective entropy 
barrier to extinction. The proportionality factor is the 
absolute value of the area between the activation trajec- 
tory and relaxation trajectory [l4| . see an example in 
Fig. [31 In contrast to scenario A, the pre-exponential 
factors in Eq. ([93]) are N- independent. Using Eqs. ([20]) . 
([2~T]) and ([9H]) . one can rewrite Eq. (|M]) in a more concise 
form: 

2tt 



aH pq (q 1 ,0) 



JV[S">( gi )-S<«(g 2 )] + [*( 9l )-*(<, 2 )] 



(94) 



As mentioned above, determining the MTE in scenario 
B does not require any information about the small-n 
recursive solution. Furthermore, Eq. (|93p formally co- 
incides with the result of Escudero and Kamenev [22l |. 
who calculated a different quantity: the mean time to 
escape from one metastable state into another. Finally, 
the same result ([93]) can be also obtained for the mean 
time to escape to an absorbing state at infinity, as in the 
particular example considered by Meerson and Sasorov 
[2l| . The reason for these coincidences is that, in all 
these systems, a constant probability current sets in be- 
yond the repelling fixed point of the rate equation. It is 
the magnitude of this current, carried by the slow WKB 
mode, rather than the exact nature of the target state 
for escape (an absorbing state at zero, infinity or an- 
other metastable state), that determines, in the leading 
and subleading orders in N, the mean escape rate from 
a metastable state. 

To conclude this section, the QSD (another main result 
of this work) is given by four overlapping asymptotes: (i) 
the recursive solution ([57]) . valid for 1 < n <C N, (ii) the 
slow-mode WKB solution ([26]). valid for ni - n » N 1 / 2 
and n 3> 1, (iii) the boundary-layer solution ([84"]) . valid 
for \n — ni | <C ri\ , and (iv) the fast-mode WKB solution 
pg]). valid for n — n\ > N 1 / 2 . 
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B. Single-step processes 

For completeness, we briefly consider the special case of 
single-step processes, where only W±\{Nq) = W±(Nq) = 
Nw±(q)+u±(q)+. . . are present. Here Eq. (|9"5)) simplifies 
considerably. Performing calculations similar to those in 
scenario A (see Sec. HV[) and using Eqs. (fT6f and (|53|) and 
the fact that w + (qi^ 2 ) = W-(qi. 2 ), one obtains the MTE 



2?re 



P 2 ( 



)dq 



aw + (q 2 )^\S"( qi )\S"(q 2 ) 



(95) 



As expected, this result coincides with the single-step re- 
sult of Ref. [22| for the mean time of a population switch 
between two metastable states. 



FIG. 5: (color online). Shown are typical zero-energy trajec- 
tories of the WKB Hamiltonian (|12p in scenario B, close to 
the bifurcation where qi — q\ <C 1. Here, the activation tra- 
jectory is a parabola, and AS given by Eq. (|100p is the area 
of the shaded region. 



C. Extinction near saddle-node bifurcation point 

Here we calculate the MTE near the characteristic 
(saddle-node) bifurcation of scenario B. At the bifurca- 
tion, the nontrivial attracting fixed point q = q 2 of the 
rate equation merges with the repelling point q = q\. 
Above but near the bifurcation point q<z — qy <C 1- As 
a result, the momentum p on the activation trajectory 
is much smaller than unity, see Fig. [5] One can always 
define the parameter N such that, at the bifurcation, 
Qi = 12 — 1- Furthermore, near the bifurcation qi = 1 — 5 
and q 2 = 1 + 8, where the exact definition of 8 -C 1 will 
appear shortly. 

Let us Taylor-expand H(q,p) from Eq. (JT2"|l in the 
vicinity of q = 1 and p = 0. As we expect p a {q) to 
be ~ (q 2 — qi) 2 , we neglect the terms of order (q — l)p 2 
and higher and arrive at the following equation for the 
zero-energy phase trajectories p = p(q) close to q = 1: 



+ ±( q -l)2 rw >;(l) + ±pr 2 w r (l) 



As S"(q) = p' a (q), we find S"(q 2 ) 



S"(qi) 



28 2 \H qqp (l, 0)\/H pp (l, 0). Furthermore, the action 
SW(q) = f q p a (Z)<% is given by 

\H qqP (l,0)\ 



-j~Y^ qi 2 ' qiq2 



whereas 

AS = S^( qi )-S^(q 2 ) = 



4|g MP (l,0)| 3 
3H PP {1,0) 



(99) 



(100) 



is the area of the shaded region in Fig. [Sj 

As in scenario A, the sub-leading WKB correction van- 
ishes near the bifurcation. Indeed, we Taylor-expand 
Eq. (1151) in the vicinity of q = 1 and p — 0, keep only 
leading-order terms, and obtain 

X>r(l)Si(g) - \r 2 w r {l){q- 1K'(1) - r<(l) ~ 0. 

r 

(101) 

0. (96) Using Eq. (IMj) . we find that the second and third terms 
cancel out, and so s[ (q) can be chosen zero. As a result, 



As can be checked a posteriori, the terms in Eq. scale Eq. yields the MTE near the bifurcation: 
as follows: H P {1,0) - <5 2 , H qp (l,0) - S 2 , H qqp (l,0) 



O(l), and H pp (l,0) = 0(1). Therefore, the term 
(q — l)H qp (l,0) ~ <5 3 can be neglected. The nontriv- 
ial solution of Eq. (|96|) yields the activation trajectory 
p = p a (q)- a parabola with the roots q = qi and q = q 2 . 
To simplify the notation, we use Eq. (I16p and evaluate 
the small difference f/2 — 9i by expanding the algebraic 
equation H p (q, 0) = in the vicinity of q = 1. Neglecting 
the term {q — l)H qp (l,0), we obtain 



2tt 



a\H qqp (l,0)\5 



exp -ND 2 B S : 



(102) 



<?2 - gi 



2g p (l,0) 
1^(1,0)1 ' 



(97) 



where H qqp (qi,0) < 0. The activation trajectory can be 
written as 

I^ OTP (1,0)| 



p a (q) = 



H pp (l,0) 



(q-qi)(q- 52) 



(98) 



where H qqp (l,0) and £> B = ^1^(1,0)1/^(1,0) = 
O(l) should be evaluated at the bifurcation. The appli- 
cability criterion of this result is iV<5 3 ^ 1. For suffi- 
ciently large N this strong inequality is compatible with 
the strong inequality S <§; 1 which describes closeness to 
the bifurcation. At S - iV -1 / 3 Eq. (fT02|) predicts the 
following scaling of the MTE with N: r ~ N 1/3 /a. 

We notice that Eq. (|102[) docs not require any infor- 
mation about the QSD in the region of small n. Indeed, 
as was mentioned in section V A, the exact nature of the 
target state is of no significance here. Note that the same 
scaling of the effective entropy barrier with the distance 
from the bifurcation appears in the context of escape 
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from one metastable state to another [l2j. Finally, the 
same scaling near the bifurcation is observed in continu- 
ous systems, driven by external delta-correlated gaussian 
noise and therefore describable by a Fokker-Planck equa- 
tion m n2 . 



D. Example: reactions 2A 3A and A — > 



Let us illustrate the extinction scenario B on the fol- 
lowing set of reactions: binary reproduction 2A — > 3A, 
the reverse process 3A A 2A, and linear decay A A 0. 
Here 

Wi = Mn-l) ^ + an(n-l)(n-2) _ 

2 ' 6 

Rescaling time pit — > t, and denoting 7 = 1 — <5 2 , <5 2 = 
l-8cr/i/(3A 2 ) > 0, and N = 3A/(2ct), we arrive at Eq. © 
with 

g 3 2q 2 
w-i = h q , 101 = — , 

7 7 

u_! = , u 1 = -^. (104) 

7 7 

In the rescaled notation the fixed points are q = 
(attracting point), gi = 1 — 8 (repelling point), and 
q2 = 1 + 5: another attracting point around which the 
metastable population resides. The WKB Hamiltonian 
(fT2"j) takes the form 



H(q,p) = (^- + q) (e- p -l) + ^-(e"- 1) 
V 7 / 7 

Solving the equation H[q,p a (q)] = yields 



Pa(q) = S (q) = In ( - + — 



Therefore, 



9 

g" -7 
<7 2 (<7 + 7) 



(105) 



(106) 



(107) 



= -5/(1 - 8), and S"(q 2 ) = 8/(1 + 8). Further- 
more, using Eq. (|104[) we obtain 



In 



'92 

and 



= 2 [ (5 — VT — d 2 arctan 



exp 



dq 



1 + 8 
1-5' 



(108) 



Plugging everything into Eq. (|95j) yields the MTE in 
physical units: 



(109) 
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FIG. 6: (color online). The natural logarithm of the ex- 
tinction rate E — 1/r versus 5 for N = 200 for the reac- 
tions A A 0, 2 A 4 3 A and 3^4, A 2 A. The analytical 
result (|109[l (the solid line) is compared with the quantity 
— [ln(l - Po um (t))]/t (the asterisks) extracted from a numer- 
ical solution of (a truncated) master equation Eq. ((2| with 
rates (fTU5]l . 



Here 

AS = 2 (5 - y/l - 5 2 arctan 6 ^ (110) 

is a monotone increasing function of 5] its asymptotes are 

f (2/3)<5 3 + (4/15)<5 5 + . . . , <5«1, 
\ 2 - 7iy2(l -<$) + . . . , 1 - d « 1 . 

Near the bifurcation, 8 <C 1, Eq. (|109p becomes 



7T / 2 , 

r = — exp -iV<5 3 



(111) 



This is in agreement with Eq. (I102p . Indeed, at the bi- 
furcation one has H pp (l,0) = 4, H qqp (l,0) — —2, and 
so D% = 1/2. Equation (fTTT| is valid when NS 3 > 1. 
This inequality, combined with <5 1, yields the double 
inequality TV -1 / 3 <5<1. 

We compared our analytical result (|109p with numer- 
ical solutions of (a truncated) master equation © with 
rates ffT03]> at N = 200 and different values of 5. The 
comparison is presented in Fig. Very good agreement 
is observed for not too small 8, when the effective entropy 
barrier NAS is sufficiently high. 

To determine the QSD in this example, one needs to 
find S (/) (<7) and s[ f) (q) from Eqs. (|TI|) and ([17]): 



S if) (q) = 2^7 arctan ( 

\V7 



In 



7 + g^ 
2g 



hi 



7 + g^ 

„5/2 



- 1 

(112) 



Then Eqs. ([24] ) . (|26 j) . (|37 ]) and (|84j ) yield the QSD in 
terms of four overlapping asymptotes. The QSD and its 
comparison with numerics are shown in Fig. and very 
good agreement is observed. 
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cesses, although for such processes our general results 
simplify considerably. The results also simplify near the 
characteristic bifurcations of scenarios A and B. A num- 
ber of previous results for the mean time to extinction 
follow from our equations in particular cases. 

We have observed that, in models belonging to extinc- 
tion scenario B, the mean time to population extinction 
formally coincides, in the leading and subleading orders 
in N , with the mean time to population escape in stochas- 
tic population models where the metastable population 
switches to another metastable state or to an absorbing 
state at infinity. The reason for this coincidence is that, 
in all these systems, a constant probability current sets in 
beyond the repelling fixed point of the rate equation. It 
is the magnitude of this current, rather than the precise 
nature of the target state for escape, which determines, in 
the leading and subleading orders in N, the mean escape 
rate from a metastable state. 

The situation is quite different in extinction scenario A. 
Here the exact nature of the target state (the absorbing 
state at n = 0) or, more precisely, the rate constants of 
the effective linear branching processes at small n, affects 
the pre-exponent of the mean time to extinction. 



FIG. 7: (color online). (a): the QSD tt„ for 8 = 1/3 

and N = 10 3 for the reactions A A 0, 2A 4 3A and 
3A — > 2A. The QSD includes four overlapping asymptotes: 
the fast-mode solution (FM), the slow-mode solution (SM), 
the boundary-layer solution (BL), and the small-n asymptote, 
(b) a comparison between the QSD in (a) (the solid line) and 
a numerical solution of (a truncated) master equation Q with 
rates (|103[) (the dashed line). 
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VI. SUMMARY 

This work dealt with extinction of an isolated long- 
lived stochastic population describable by a continuous- 
time Markov process. We have identified two generic ex- 
tinction scenarios, A and B, based on the stability prop- 
erties of the fixed points of the population size dynam- 
ics, predicted by the rate equation. For each of the two 
scenarios we have calculated the mean time to extinc- 
tion (MTE) and the quasi-stationary probability distri- 
bution (QSD). The calculations involve a systematic use 
of WKB method, where 1/N, the typical inverse popu- 
lation size in the metastable state, serves as a small pa- 
rameter of the theory. The WKB theory is supplemented 
by two additional approximations in the regions where it 
breaks down. One of them is a small-n expansion of the 
quasi-stationary master equation which brings about a 
recursion relation (in both scenarios A and B). The sec- 
ond one is the Fokker-Planck equation, obtained via the 
van Kampen system-size expansion. The latter is valid 
in small regions around the non-trivial fixed points of the 
rate equation (in scenario B). 

The theory is not limited to single-step stochastic pro- 



Appendix A 

Here we derive Eq. (|55|) for the arbitrary constants C; 
appearing in Eq. For that we need to solve the 

following set of K + 1 linear equations: 

K 

Co + ^2K mc i = fm, m=l,2,...,K + l, (Al) 

i=l 

where the constants f%, fy, . . . , /k+i can be expressed via 
/i by using the recursive Eq. P0|) with n — 1,2, ... ,K 
(assuming that fj = for j < 1). 

As before, we denote the roots of Eq. (J3TJ) by 
Ao, Ai, . . . , \k, where Ao = 1, and define the following 
quantities: 



s\ ' = Ai + A 2 + A 3 + • • • + Xk 

s 2°' — ^1^2 + A2A3 + A3A4 + • • • + Xk-iX-K 
sf^ — A1A2A3 + A1A3A4 + A1A2A4 + • ■ ■ 



s%> = A1A2-.-A.pf. (A2) 
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The superscript (0) means that the K roots which con- 
tribute to s\ are all of the roots of Eq. (j3"Tj) except Xq. 
In the same manner we can define for 1 < i < K 



Let us apply this formula to Eq. which has exactly 
the roots Ai,A2,--- ,Xk- For convenience, we rewrite 
Eq. d33J as 



(*) 

(0 
»i 



(0 



1 

1 + Ai + A 2 + • • • + Ai_i + A J+1 + • • • + X K 

K 

AjA m 

< j < TO < if 



— A1A2 • • • Ai_iA,;+i • • • Xk , 



(A3) 



where the root Ai does not contribute. 

To obtain the coefficient Co we multiply the first of 
Eqs. (|Alj) by the second by — and finally the 
last by (—l) K Sft\ By adding the equations we obtain 



f «(°) f 4_ a (°) f 



+ (-1)^/a + i 



= C 
+Ci 



1 - s 



(0) , JO) 



-... + (-!)*« 



A" 



A r l -ArM 0) + Ar 3 4 0) ---- + (-i) A 'Ar A "- 1 ^ ) 



c 



A 



A A _ A A s i H r(. _J -J A A S A 



(A4) 



By using Eqs. (|A2j) we can rewrite the coefficient of Co 
in the right hand side of Eq. (|A4[) as 



l_ a (o) +fl (o)_. . .+(_i)« s ^ = (i_ Al )(l-A 2 ) • • • (1-X K ). 

(A5) 

Furthermore, the coefficient of C in Eq. (|A4j) satisfies 



A^-A^ + A- 3 ^ 



+ (-l)"A 



(Ai - Ai)(Ai - A 2 ) •• • (Ai - A,) •• • (Ai - X K ) 



X 



K+l 



(A6) 



Clearly, this expression vanishes for alH > 1. Therefore, 
Eqs. dSU and (JMJ yield 



A" 



A 



1-A^<(0)-A 2 ^ W ;(0) A A '«4(0) =0. (A9) 

r=l r=2 

In terms of the coefficients a, we have ao = 1, ai = 

-K(o) + ---+^(o)],a 2 = -K(o) + ---+^(o)],---, 

and a K = -w' K (Q). Therefore, using Eqs. (jA"7|) and 
we can rewrite Co as 

^ _ /a+i - feKj + ■ ■ ■ + w' K (Q)] f lW ' K (0) 

°° - A ' 



1=1 



(A10) 



Now, using Eq. (|30|) with n = K, we obtain a relation 
between and fj<x- Plugging it into (|A10[) we have 



Co 



{f K -f K -M(o)+---+<m 

^na-Ai) 



i=l 



/iK-i(fl)+4(o)]} 



(All) 



One can use this argument repeatedly if — 1 more times, 
and obtain that the expression in the curly brackets in 
Eq. (|Allj) equals f x . In addition, by virtue of Eq. (|A"8)| 
w' K (0) satisfies 



y A-(°) = ~ a K = 



AiA 2 • • • Ajf 



(A12) 



Therefore, Co is given by 

A 



(-1)^/1 II A * 



C n = 



A 



A 



(A13) 



na-Ai) n^-v 1 ; 



i=l 



f B (°) f _i_ e (°) f _i_ f „(0) f 

Co = Jl- S l J2 + S 2 h b (, — lj S K JK+1 ^ A7 ^ 



na-A.) 



i=l 



To calculate the numerator of Eq. (|A7J) we use Viete's 
formula. Given a polynomial equation 



clkX + ■ ■ ■ + a\x + ao = , 

whose roots are Ai,A2,--- , A^, the expressions sf\ 
given by Eq. (|A2|> . satisfy 



J0) 



•1) 



a K 



(A8) 



thereby proving Eq. (|35|) for i = 0. 

This proof can be generalized to the rest of the coeffi- 
cients Cj. Here one has to multiply the first of Eqs. (jAlj) 
by Sq , the second by — s± and finally the last by 
( — l) K sS, and add all the equations. This yields 



71 - s l J2 + S 2 J3 



= C 
+Cr 



(0 , „(0 



A„« 



-...+(-!)*« 



A 



Ar^Ar^ + Ar 3 ^ 



c 



A 



x^-x-^ + ... + (-i) K x- K -^ 



(A14) 
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The coefficient of C, in the right hand side of Eq. (|A14|) 
satisfies 



A 7 



K^ + K' 6 ^ -••• + (-!)* A 



A" 



(Ai-l)(Aj-Ai) ■ ■ ■ (Ai-Aj_i)(A«-A i+ i) ■ ■ ■ (Aj-A g ) 
Af+1 

(A15) 



because A,; is absent from Sj, see Eq. (IA3D . On the 



other hand, the coefficients of all other Cj^i in the right 
hand side of Eq. (|A14l) can be shown to be equal zero. 
Therefore, 



A 



K+l 



Ci 



h-sfh + sfh 



n (^- A i) 

i = o 

3 + i 

(A16) 

By using the recursive equation ([50)) repeatedly and by 
using Eq. (|A3|) . one finally obtains Eq. ((331) . 

Finally, Co from Eq. (IA13|) can be expressed through 
the reaction rate constants w[(0) , w' 2 (0) , ■ ■ ■ ,w' K (0), see 
Eq. ((33)) . Indeed, let us expand the denominator of 
Eq. fAHfl) and use Eq. (|Al|) : 



(i - Ar x )(i - a^ 1 ) • • • (i - a 7 /) = i - E v 1 



1 > J 



Ai • • • Xjc 



(-l)^ + (-l) 



K-l.(O) 
A if-1 



f 1 . 

(A17) 



Using Eqs. (IA"8]) and (|A12|) . we rewrite Eq. (|A17|) as 

(i - A r x )(i - k 1 ) ■■■(i- A^ 1 ) = 1 - K(0) + • • • 

+w' K (0)] - K(0) + • • • + w' K (0)] <(0) 

= l-w[ (0) - 2w 2 (0) - 3v/ 3 (0) (0) . ( A18) 

Plugging this into Eq. (jA13p one obtains Eq. ((3d]) . 



Appendix B 

Here we show that, in extinction scenario A, the two 
real and positive roots of Eq. (j3"Tj) are Ao = 1 and 

< Ai < 1, whereas all other roots obey the inequal- 
ity |Ai>i| > Ai. We start by showing that the posi- 
tive root of Eq. (133)) obeys the inequalities < A x < 
1. The left hand side of Eq. ((35)) side is a monotone 
decreasing function of A. At A = 1 it is equal to 

1 - w[(0) -2w' 2 (0) Kw' K (0). This quantity is 

negative, as n = is a repelling fixed point, so the 
rescaled reaction rate constants satisfy the inequality 



w[(0) + 2w' 2 (0) + ■ ■ ■ + Kw' K (0) - 1 > 0. On the other 
hand, at A = the left hand side is 1. Hence, < Ai < 1. 

Now we will prove by contradiction that all other (neg- 
ative or complex) roots satisfy the inequaltity |Aj>i| > 
Ai. Assume by contradiction that there exists a root Xj 

Then by assump- 



ae 



i0 



so that |Aj| < Ai. Denote Ac- 
tion a < Ai. Substituting Xj into Eq. ((33J we have 



1-«j1(0)A 3 -- 
= 1 — vu'i (0)a cos ( 
-w' K (0){acos9 + 



K 



(0)(A j 



+ Af ) 



a K cos K6) + i(- 



0,(B1) 



where both real and imaginary parts have to vanish sep- 
arately. Now we substitute Ai into Eq. ((33)) and use 
Eq. djjfll 



= 1 - 10^(0) Ai - 

< l-w[(0)a--- 

< 1 — w' 1 (0)a cos ( 
- w' K (0)(acos6 4 



-w' K (0)(\i + 



K 



(0)(a- 



+ Af) 



+ a K ) 



,K 



cosK9)=0, (B2) 



where the last inequality holds since Xj is complex or 
negative, so 9 =/= 0, and there exists some m for which 
cos m9 < 1. 

Eq. ()B2j) shows a contradiction < 0, so all roots obey 
|Aj>i| > Ai. As a result, at n > 1, the recursive solution 
(157)1 reduces to Eq. (f55|). where A 1 is given by Eq. (1551) . 

Finally, we show that A\ > 0, and so the asymp- 
tote ((55)) is always positive. First, by using Eq. (|A12[) . 
one can see that the numerator in Eq. ((55)) is always 
negative. What is the sign of the denominator? Here 
Ai — A < 0, whereas all other terms in the product are 
positive. Indeed, for Xj < one has Ai — Xj > 0. For 
any complex A^, there is also a complex conjugate root 
Xk = Xj. Therefore, by writing Xj = a + ib, one has 
(Ai -Xj)(Xi - Xj) = (Ai -a) 2 + b 2 > 0. So, A 1 is always 
positive, and so is the asymptote 7r(g) from Eq. 



Appendix C 

Here we show that, in extinction scenario B, the root 
Ao = 1 of Eq. ((5T)) has the smallest absolute value among 
all of the roots. To this end we will prove that the roots 
of Eq. ((55)) obey the inequality |Aj>q| > 1. Let us denote 
Aj>o = a,je l9j , in general a complex number, and assume 



by contradiction that a,j 
Eq. ((33)) . we obtain 



< 1. Plugging A = Xj into 



1 - aj cos(6 j )w' 1 (0) - [a,j cos(9j)+a 2 cos(29j)]w' 2 (0) 



a,- cosl 



) + ••• + af cos(K6j)]w' K (0) + i (•••) = 0, 



(CI) 



where the real and imaginary parts must vanish sep- 
arately. As we have assumed cij < 1, we have 
aj 1 cos(m9j) < 1 for all integer m > 0. Therefore, we 
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can write for the real part of Eq. (|C1|) : 

= 1 - a,j cos(0j)ii>i(O) - [ aj cos(6»j) + aj cos{26 1 )]w' 2 (0) 

[aj cos(9j) + ■ ■ ■ + af cos(K6j)]w' K (0) 

> l-w[(0)- 2w' 2 (0) Kw' K (0) > , (C2) 



where the last inequality follows from n — being an 
attracting fixed point of the rate equation. Equation (|C2j) 
shows a contradiction > 0. Hence aj > 1, and all the 
roots of Eq. ([55)) obey the inequality |A,>q| > 1. 



M.S. Bartlett, Stochastic Population Models in Ecology [17 
and Epidemiology (Wiley, New York, 1961). 
S. R. Beissinger and D. R. McCullough (Editors), Pop- [18 
ulation Viability Analysis (University of Chicago Press, 
Chicago, 2002). [19 
H. Andersson and T. Britton, Stochastic Epidemic Mod- 
els and Their Statistical Analysis, Lect. Notes Stat., Vol. [20 
151 (Springer, New York, 2000). 

M.S. Samoilov and A. P. Arkin, Nature Biotech. 24, 1235 [21 
(2006); M. Assaf and B. Meerson, Phys. Rev. Lett. 100, 
058105 (2008). [22 
N.G. van Kampen, Stochastic Processes in Physics and 
Chemistry (North- Holland, Amsterdam, 2001). [23 

C. W. Gardiner, Handbook of Stochastic Methods 
(Springer Verlag, Berlin, 2004). 

We will assume that the stochastic population does not [24 
exhibit an unlimite d g rowth (escape to infinite popula- [25 
tion size), see Ref. [2l[, and n = is the only absorbing 
state. [26 
M. Assaf and B. Meerson, Phys. Rev. E 74, 041115 
(2006). 

M. Assaf and B. Meerson, Phys. Rev. Lett. 97, 200602 
(2006). [27 
CM. Bender and S.A. Orszag, Advanced Mathematical 
Methods for Scientists and Engineers (Springer, New [28 
York, 1999). [29 
R. Kubo, K. Matsuo, and K. Kitahara, J. Stat. Phys. 9, [30 
51 (1973). 

M.I. Dykman, E. Mori, J. Ross, and P.M. Hunt, J. Chem. 
Phys. 100, 5735 (1994). 

V. Elgart and A. Kamenev, Phys. Rev. E 70, 041106 [31 
(2004). 

V. Elgart and A. Kamenev, Phys. Rev. E 74, 041101 [32 

(2006) . 

M. Assaf, A. Kamenev, and B. Meerson, Phys. Rev. E 
78, 041123 (2008). 

D. A. Kessler and N.M. Shnerb, J. Stat. Phys. 127, 861 

(2007) . 



B. Gaveau, M. Moreau, and J. Toth, Lett. Math. Phys. 
37, 285 (1996). 

C. R. Doering, K.V. Sargsyan, and L.M. Sander, Multi- 
scale Model, and Simul. 3, 283 (2005). 

M. Assaf and B. Meerson, Phys. Rev. E 75, 031122 
(2007). 

J.W. Turner and M. Malek-Mansour, Physica A 93, 517 
(1978). 

B. Meerson and P.V. Sasorov, Phys. Rev. E 78, 
060103(R) (2008). 

C. Escudero and A. Kamenev, Phys. Rev. E 79, 041149 
(2009). 

At nontrivial fixed points of the rate equation, q = qi = 
m/N > 0, the two real roots of the equation H(q,p) = 
merge at p = 0. 

If S"'(0) = 0, then it is more convenient to use Eq. (|17|1 . 
J.N. Darroch and E. Seneta, J. Appl. Probab. 4, 192 
(1967). 

We remind the reader that, in order to obtain Eq. ©, 
we rescaled the reaction rates and time by the linear de- 
cay rate constant w'_i(0) = a. To express the MTE in 
physical units one needs to put the factor a back. 
M. Assaf, A. Kamenev, and B. Meerson, Phys. Rev. E 
79, 011127 (2009). 

I. Nasell, J. Theor. Biol. 211, 11 (2001). 

O. Ovaskainen, J. Appl. Prob. 38, 898 (2001). 

In Refs. p8l. |2^] the parameter Ro was defined as A (TV — 

l)//i. Therefore, their result for the MTE looks slightly 

different, but it is actually identical to Eq. (|71|) in the 

leading and subleading orders of N. 

M.I. Dykman and M.A. Krivoglaz, Physica A 104, 480 
(1980). 

The distance to bifurcation 6i is defined in Refs. [12| and 
[3~D ] differently than in the present work. Their parameter 
Si is related to our parameter S as 5i = 5 2 , so the entropy 
barriers in Refs. [l^] and [3l[ scale with Si as &\^ 2 . 



